library(ggplot2)
library(GGally)
library(wesanderson)
library(mgcv)
library(tidymv)
library(nnet) # neural networks/multinomial logistic regressions
library(AER) ## coeftest from nnet
library(readxl)
library(stargazer)
library(MNLpred)
library(survival)
library(survminer)
library(contsurvplot)
library(modelsummary)
library(marginaleffects)
library(patchwork)
#
## Figure A1. PSDI with density ####
edi2d <- ggplot(parsys_vdem, 
                aes(x = v2ps_democracy_adj, y = v2x_polyarchy)) + 
  geom_point()+
  stat_density_2d(aes(fill = after_stat(density)), geom = "raster", contour = FALSE) +
  scale_fill_continuous(type = "viridis") +
  xlab("Party-System Democracy Index")+
  ylab("Electoral Democracy Index (V-Dem, 2023)") +
  labs(title = "V-Dem")+
  geom_hline(yintercept = 0.5, color = "black", linetype = "dashed") +
  scale_y_continuous(expand = c(0,0), breaks=seq(0, 1, by = .1)) +
  scale_x_continuous(expand = c(0,0), breaks=seq(-1, 1, by = 0.1)) +
  stat_cor(label.y = .05, label.x = 0.15, color ="white", r.digits = 3, p.accuracy = 0.001)+
  theme_light()+
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 

p52d <- ggplot(parsys_p5, 
               aes(x = v2ps_democracy_adj, y = polity2)) + 
  geom_point()+
  stat_density_2d(aes(fill = after_stat(density)), geom = "raster", contour = FALSE) +
  scale_fill_continuous(type = "viridis") +
  xlab("Party-System Democracy Index")+
  ylab("Polity Score (Polity5, 2018)") +
  labs(title = "Polity5")+
  geom_hline(yintercept = 0, color = "black", linetype = "dashed") +
  scale_y_continuous(expand = c(0,0), breaks=seq(-10, 10, by = 1)) +
  scale_x_continuous(expand = c(0,0), breaks=seq(-1, 1, by = 0.1)) +
  stat_cor(label.y = -9.5, label.x = 0.15, color ="white", r.digits = 3, p.accuracy = 0.001)+
  theme_light()+
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fh2d <- ggplot(parsys_fh, 
               aes(x = v2ps_democracy_adj, y = Total)) + 
  geom_point()+
  stat_density_2d(aes(fill = after_stat(density)), geom = "raster", contour = FALSE) +
  scale_fill_continuous(type = "viridis") +
  xlab("Party-System Democracy Index")+
  ylab("Freedom in the World (FreedomHouse, 2023)") +
  geom_hline(yintercept = 50, color = "black", linetype = "dashed") +
  scale_y_continuous(expand = c(0,0), breaks=seq(0, 100, by = 10)) +
  scale_x_continuous(expand = c(0,0), breaks=seq(-1, 1, by = 0.1)) +
  labs(title = "Freedom House")+
  stat_cor(label.y = 5, label.x = 0.15, color ="white", r.digits = 3, p.accuracy = 0.001)+
  theme_light()+
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 

edi2d + p52d + fh2d

ggsave(filename = here("Figures/PSDI_democracy2d.pdf"),
       dpi = "retina",
       width = 14,
       height = 8)

## Figure A2. Western Europe and North America ####
wena_psdi <- ggplot(subset(parsys_vdem, e_regionpol_6C ==5),
                   aes(year)) +
  geom_line(aes(y = v2ps_democracy_adj, color = "PSDI"))+
  geom_line(aes(y = v2ps_democracy_codehigh_adj, color = "PSDI"), linetype="twodash")+
  geom_line(aes(y = v2ps_democracy_codelow_adj, color = "PSDI"), linetype="twodash")+
  xlab("Party-System Democracy Index")+ 
  ylab("Election Year") +
  scale_y_continuous(breaks=seq(0, 1, by = 0.25)) +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  theme_light()+
  labs(title = "Western Europe and North America")+
  theme(legend.position = "bottom")+
  facet_wrap(~country_name, scale="fixed")

## Figure A3. Eastern Europe and Central Asia ####
eeca_psdi <- ggplot(subset(parsys_vdem, e_regionpol_6C.x ==1),
                   aes(year)) +
  geom_line(aes(y = v2ps_democracy_adj, color = "PSDI"))+
  geom_line(aes(y = v2ps_democracy_codehigh_adj, color = "PSDI"), linetype="twodash")+
  geom_line(aes(y = v2ps_democracy_codelow_adj, color = "PSDI"), linetype="twodash")+
  xlab("Party-System Democracy Index")+ 
  ylab("Election Year") +
  scale_y_continuous(breaks=seq(0, 1, by = 0.25)) +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  theme_light()+
  labs(title = "Eastern Europe and Central Asia")+
  theme(legend.position = "bottom")+
  facet_wrap(~country_name, scale="fixed")

## Figure A4. Latin America and the Caribbean #### 
lac_psdi <- ggplot(subset(psdi.data.b, e_regionpol_6C ==2),
                  aes(year)) +
  geom_line(aes(y = v2ps_democracy_adj, color = "PSDI"))+
  geom_line(aes(y = v2ps_democracy_codehigh_adj, color = "PSDI"), linetype="twodash")+
  geom_line(aes(y = v2ps_democracy_codelow_adj, color = "PSDI"), linetype="twodash")+
  xlab("Party-System Democracy Index")+ 
  ylab("Election Year") +
  scale_y_continuous(breaks=seq(0, 1, by = 0.25)) +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  theme_light()+
  labs(title = "Latin America and Caraibbean")+
  theme(legend.position = "bottom")+
  facet_wrap(~country_name, scale="fixed")

## Figure A5. Middle East and North Africa ####

parsys_vdem_mena <- psdi.data.b%>%
  filter(country_name != "Oman" &  
           country_name != "United Arab Emirates")

mena_psdi <- ggplot(subset(parsys_vdem_mena, e_regionpol_6C==3),
                   aes(year)) +
  geom_line(aes(y = v2ps_democracy_adj, color = "PSDI"))+
  geom_line(aes(y = v2ps_democracy_codehigh_adj, color = "PSDI"), linetype="twodash")+
  geom_line(aes(y = v2ps_democracy_codelow_adj, color = "PSDI"), linetype="twodash")+
  xlab("Party-System Democracy Index")+ 
  ylab("Election Year") +
  scale_y_continuous(breaks=seq(0, 1, by = 0.25)) +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  theme_light()+
  labs(title = "Middle East and Northen Africa")+
  theme(legend.position = "bottom")+
  facet_wrap(~country_name, scale="fixed")

## Figure A6. Sub-Saharan Africa ####
ssa_psdi <- ggplot(subset(parsys_vdem, e_regionpol_6C ==4),
                  aes(year)) +
  geom_line(aes(y = v2ps_democracy_adj, color = "PSDI"))+
  geom_line(aes(y = v2ps_democracy_codehigh_adj, color = "PSDI"), linetype="twodash")+
  geom_line(aes(y = v2ps_democracy_codelow_adj, color = "PSDI"), linetype="twodash")+
  xlab("Party-System Democracy Index")+ 
  ylab("Election Year") +
  scale_y_continuous(breaks=seq(0, 1, by = 0.25)) +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  theme_light()+
  labs(title = "Sub-Saharan Africa")+
  theme(legend.position = "bottom")+
  facet_wrap(~country_name, scale="fixed")

## Figure A7. Asia Pacific ####
ap_psdi <- ggplot(subset(parsys_vdem, e_regionpol_6C ==6),
                 aes(year)) +
  geom_line(aes(y = v2ps_democracy_adj, color = "PSDI"))+
  geom_line(aes(y = v2ps_democracy_codehigh_adj, color = "PSDI"), linetype="twodash")+
  geom_line(aes(y = v2ps_democracy_codelow_adj, color = "PSDI"), linetype="twodash")+
  xlab("Party-System Democracy Index")+ 
  ylab("Election Year") +
  scale_y_continuous(breaks=seq(0, 1, by = 0.25)) +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  theme_light()+
  labs(title = "Asia Pacific")+
  theme(legend.position = "bottom")+
  facet_wrap(~country_name, scale="fixed")

## Figure A8. PSDI and RoW ####
row.plot <- ggplot(subset(parsys_vdem, !is.na(v2x_regime)), 
         aes(x = v2ps_democracy_adj, y = factor(v2x_regime), 
             fill = factor(v2x_regime)))+ 
  geom_density_ridges(
    jittered_points = TRUE, quantile_lines = TRUE, scale = 0.9, alpha = 0.7,
    vline_size = 1, vline_color = "red",
    point_size = 0.4, point_alpha = 1,
    position = position_raincloud(adjust_vlines = TRUE))+
  scale_fill_manual(values = wes_palette("BottleRocket2", n = 4),
                    name = "Regimes of the World",labels = c("Closed Autocracies", "Electoral Autocracies",
                                                             "Electoral Democracies", "Liberal Democracies"))+
  xlab("Party-System Democracy Index")+
  ylab("Regimes of the World (V-Dem, 2023)") +
  theme_light()+
  theme(axis.text.y =element_blank(), 
        legend.position = "bottom") 

## Figure A9. PSDI Disaggregating 4 cases ####
cases_app <- parsys_vdem %>%
  filter(country_name == "Germany"| country_name =="Mexico"
         |country_name == "India" |country_name == "Russia")

ld.psdi_app <- ggplot(cases_app,
                      aes(year, 
                          group = factor(country_text_id.x), 
                          color = factor(country_text_id.x))) +
  geom_line(aes(y = 1-v2ps_democracy_government), color = "steelblue",  linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_government_codehigh), color = "steelblue", linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_government_codelow), color = "steelblue", linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_adj), color = "darkgreen",  linetype="twodash")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_codehigh_adj), color = "darkgreen", linetype="twodash")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_codelow_adj), color = "darkgreen", linetype="twodash")+
  geom_line(aes(y = v2x_polyarchy), color = "firebrick", linetype="solid")+
  gghighlight(cases_app$country_name =="Germany", 
              unhighlighted_colour = alpha("lightblue", 0),
              use_direct_label = F,
              label_key = country_name,
              label_params = list(size = 5))+ 
  xlab("Election Year")+ 
  ylab("Party-System Democracy Index") +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  labs(title = "Germany")+
  theme_light()+
  theme(legend.position = "bottom")

aut.psdi_app <-  ggplot(cases_app,
                        aes(year, 
                            group = factor(country_text_id.x), 
                            color = factor(country_text_id.x))) +
  geom_line(aes(y = 1-v2ps_democracy_government), color = "steelblue",  linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_government_codehigh), color = "steelblue", linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_government_codelow), color = "steelblue", linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_adj), color = "darkgreen",  linetype="twodash")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_codehigh_adj), color = "darkgreen", linetype="twodash")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_codelow_adj), color = "darkgreen", linetype="twodash")+
  geom_line(aes(y = v2x_polyarchy), color = "firebrick", linetype="solid")+
  gghighlight(cases_app$country_name =="India", 
              unhighlighted_colour = alpha("lightblue", 0),
              use_direct_label = F,
              label_key = country_name,
              label_params = list(size = 5))+ 
  xlab("Election Year")+ 
  ylab("Party-System Democracy Index") +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  labs(title = "India")+
  theme_light()+
  theme(legend.position = "bottom")

demo.psdi_app <-  ggplot(cases_app,
                         aes(year, 
                             group = factor(country_text_id.x), 
                             color = factor(country_text_id.x))) +
  geom_line(aes(y = 1-v2ps_democracy_government, color ="Government"),  linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_government_codehigh), color = "steelblue", linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_government_codelow), color = "steelblue", linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_adj, color = "Opposition"),  linetype="twodash")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_codehigh_adj), color = "darkgreen", linetype="twodash")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_codelow_adj), color = "darkgreen", linetype="twodash")+
  geom_line(aes(y = v2x_polyarchy), color = "firebrick", linetype="solid")+
  scale_color_manual(name = "Legend PSDI", values = c("Government" = "steelblue","Opposition" = "darkgreen"))+
  gghighlight(cases_app$country_name =="Mexico", 
              unhighlighted_colour = alpha("lightblue", 0),
              use_direct_label = F,
              label_key = country_name,
              label_params = list(size = 5))+ 
  xlab("Election Year")+ 
  ylab("Party-System Democracy Index") +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  labs(title = "Mexico")+
  theme_light()+
  theme(legend.position = "bottom")

ea.psdi_app <-  ggplot(cases_app,
                       aes(year, 
                           group = factor(country_text_id.x), 
                           color = factor(country_text_id.x))) +
  geom_line(aes(y = 1-v2ps_democracy_government), color = "steelblue",  linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_government_codehigh), color = "steelblue", linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_government_codelow), color = "steelblue", linetype="dashed")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_adj), color = "darkgreen",  linetype="twodash")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_codehigh_adj), color = "darkgreen", linetype="twodash")+
  geom_line(aes(y = 1-v2ps_democracy_opposition_codelow_adj), color = "darkgreen", linetype="twodash")+
  geom_line(aes(y = v2x_polyarchy, color ="EDI"), linetype="solid")+
  scale_color_manual(name = "", values =c("EDI" = "firebrick"))+
  gghighlight(cases_app$country_name =="Russia", 
              unhighlighted_colour = alpha("lightblue", 0),
              use_direct_label = F,
              label_key = country_name,
              label_params = list(size = 5))+ 
  xlab("Election Year")+ 
  ylab("Party-System Democracy Index") +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  labs(title = "Russia and former Soviet Union")+
  theme_light()+
  theme(legend.position = "bottom")

(ld.psdi_app+aut.psdi_app)/(demo.psdi_app+ea.psdi_app)

## Figure A10. Correlation PSDI, Gvt, Opp ####
small_parsys_vdem <- parsys_vdem%>%
  dplyr::select(v2ps_democracy_adj, v2ps_democracy_government_inverted, 
                v2ps_democracy_opposition_adj_inverted,e_regionpol_6C.x)

ggpairs(subset(small_parsys_vdem, v2ps_democracy_adj>0),         # Data frame
        columns = 1:3,
        columnLabels = c("PSDI","PSDI Government","PSDI Opposition"),
        size = 0.5)+
  theme_light()

## Figure A11. Seat Share vs. Vote Share ####
vall$pavote.1 <- vall$v2pavote*0.01

psdi_vote <- vall %>%   
  dplyr::select(country_id, country_text_id, e_regionpol_6C.x,
                                      country_name, year,v2pagovsup.dum, 
                                      v2x_polyarchy.x,v2x_polyarchy_codelow.x,v2x_polyarchy_codehigh.x,
                                      pavote.1,
                                      v2xpa_antiplural) %>%
  dplyr::group_by(country_id, country_text_id, e_regionpol_6C.x,
                  country_name, year,v2pagovsup.dum,
                  v2x_polyarchy.x,v2x_polyarchy_codelow.x,v2x_polyarchy_codehigh.x) %>% 
  dplyr::summarise(
    wm_vs_autpol = weighted.mean(v2xpa_antiplural,pavote.1),
    .groups = 'drop') %>%
  as.data.frame()

psdi_vote_data <- psdi_vote %>%
  spread(v2pagovsup.dum,wm_vs_autpol) 
psdi_vote_data$v2ps_democracy_opposition_vs <- psdi_vote_data$`0`
psdi_vote_data$v2ps_democracy_opposition_adj_vs <- ifelse(is.na(psdi_vote_data$v2ps_democracy_opposition_vs), 0, 
                                                          psdi_vote_data$v2ps_democracy_opposition_vs)
psdi_vote_data$v2ps_democracy_government_vs <- psdi_vote_data$`1`
psdi_vote_data$psdi_vs <- 1-(psdi_vote_data$v2ps_democracy_government_vs+psdi_vote_data$v2ps_democracy_opposition_adj_vs)
psdi_vote_data$v2ps_democracy_vs <- rescale(psdi_vote_data$psdi_vs, to=c(0,1))
psdi_vote_data  <- psdi_vote_data%>%
  dplyr::mutate(v2ps_democracy_adj_vs = ifelse(v2x_polyarchy.x>.5, v2ps_democracy_vs, "NA")) %>%
  dplyr::mutate(v2ps_democracy_adj_vs = ifelse(v2x_polyarchy.x<=.5 & psdi_vote_data$v2ps_democracy_opposition_adj_vs ==0, 
                                               0, v2ps_democracy_vs))%>%
  as.data.frame()
psdi_vs_ss <- left_join(parsys_vdem, psdi_vote_data, by = c("country_name", "year"))

ggplot(psdi_vs_ss, 
       aes(x = v2ps_democracy_adj, y = v2ps_democracy_adj_vs, label = country_name, 
           color = factor(v2x_regime))) + 
  geom_point() +
  xlab("PSDI weighted by Seat Share")+
  ylab("PSDI weighted by Vote Share") +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, color = "black", linetype = "dashed") +
  scale_y_continuous(breaks=seq(-1, 1, by = .1)) +
  scale_x_continuous(breaks=seq(-1, 1, by = 0.1)) +
  scale_color_manual(values = wes_palette("BottleRocket2", n = 4),
                     name = "Regimes of the World",labels = c("Closed Autocracies", "Electoral Autocracies",
                                                              "Electoral Democracies", "Liberal Democracies"))+
  geom_text_repel(
    data          = subset(psdi_vs_ss, v2ps_democracy_adj ==0 & v2ps_democracy_adj_vs>0),
    max.overlaps = nrow(psdi_vs_ss),
    nudge_x      = 0.5,
    direction    = "y",
    hjust        = 0.9,
    segment.size = 0.2
  ) +
  geom_label_repel(
    data          = subset(psdi_vs_ss, v2ps_democracy_adj ==0 & v2ps_democracy_adj_vs>0),
    aes(label = year), color = "black",
    box.padding = .04,
    point.padding = .5,
    nudge_x      = 0.1,
    segment.color = "black")+
  labs(color = "RoW") +
  theme_light()+
  theme(legend.position = "bottom") 

## Table A1. List PSDI > PSDI by vote share ####
list_psdi_vs <- psdi_vs_ss%>%
  filter(v2ps_democracy_adj_vs ==0 & v2ps_democracy_adj>0)%>%
  dplyr::select(country_name, year, v2ps_democracy_adj)

xtable::xtable(list_psdi_vs, caption = "List of Observations Where PSDI is Greater than PSDI Weighted by Vote Share")

## Figure A12. PSDI and Electoral Systems ####
elsys.plot <- 
  ggplot(subset(parsys_vdem, !is.na(v2elparlel)), 
         aes(x = v2ps_democracy_adj, y = factor(v2elparlel), fill = factor(v2elparlel)))+ 
  geom_density_ridges(
    jittered_points = TRUE, quantile_lines = TRUE, scale = 0.9, alpha = 0.7,
    vline_size = 1, vline_color = "red",
    point_size = 0.4, point_alpha = 1,
    position = position_raincloud(adjust_vlines = TRUE))+
  scale_fill_manual(values = wes_palette("BottleRocket2", n = 4),
                    name = "Electoral Systems",labels = c("Others", "Mixed",
                                                          "Proportional", "Majoritarian"))+
  xlab("Party-System Democracy Index")+
  ylab("Electoral Systems (V-Dem, 2023)") +
  theme_light()+
  theme(axis.text.y =element_blank(), 
        legend.position = "bottom") 

## Figure A13. Additional Construct Validity #### 

### CHES 1999-2019 ###
CHES9919 <- read_csv("1999-2019_CHES_dataset_means(v3).csv")

# 1 = governemnt, 0 = opposition
CHES9919$govt1 <- NA
CHES9919$govt1[CHES9919$govt == 0] <- 0
CHES9919$govt1[CHES9919$govt == 0.5] <- 1
CHES9919$govt1[CHES9919$govt == 1] <- 1

# create new column with country names
CHES9919$country_name <- NA
CHES9919$country_name[CHES9919$country == 1] <- "Belgium"
CHES9919$country_name[CHES9919$country == 2] <- "Denmark"
CHES9919$country_name[CHES9919$country == 3] <- "Germany"
CHES9919$country_name[CHES9919$country == 4] <- "Greece"
CHES9919$country_name[CHES9919$country == 5] <- "Spain"
CHES9919$country_name[CHES9919$country == 6] <- "France"
CHES9919$country_name[CHES9919$country == 7] <- "Irland"
CHES9919$country_name[CHES9919$country == 8] <- "Italy"
CHES9919$country_name[CHES9919$country == 10] <- "Netherlands"
CHES9919$country_name[CHES9919$country == 11] <- "United Kingdom"
CHES9919$country_name[CHES9919$country == 12] <- "Portugal"
CHES9919$country_name[CHES9919$country == 13] <- "Austria"
CHES9919$country_name[CHES9919$country == 14] <- "Finland"
CHES9919$country_name[CHES9919$country == 16] <- "Sweden"
CHES9919$country_name[CHES9919$country == 20] <- "Bulgaria"
CHES9919$country_name[CHES9919$country == 21] <- "Czech Republic"
CHES9919$country_name[CHES9919$country == 22] <- "Estonia"
CHES9919$country_name[CHES9919$country == 23] <- "Hungary"
CHES9919$country_name[CHES9919$country == 24] <- "Latvia"
CHES9919$country_name[CHES9919$country == 25] <- "Lithuania"
CHES9919$country_name[CHES9919$country == 26] <- "Polonia"
CHES9919$country_name[CHES9919$country == 27] <- "Romania"
CHES9919$country_name[CHES9919$country == 28] <- "Slovakia"
CHES9919$country_name[CHES9919$country == 29] <- "Slovenia"
CHES9919$country_name[CHES9919$country == 31] <- "Croatia"
CHES9919$country_name[CHES9919$country == 37] <- "Malta"
CHES9919$country_name[CHES9919$country == 38] <- "Luxembourg"
CHES9919$country_name[CHES9919$country == 40] <- "Cyprus"

# create weight on left-right scale
ches.wg <- CHES9919 %>%
  group_by(country_name, year, govt1) %>%
  summarise(wm_lrgen = weighted.mean(lrgen, seat)) %>%
  as.data.frame()

# we spread the df creating one column with the weighted left-right mean
# by government (1) and opposition (0)
ches.data.wg <- ches.wg %>% 
  spread(govt1, wm_lrgen) %>%
  as.data.frame()
# create the LR index as the net difference between opposition and government
ches.data.wg$lrps <- ches.data.wg$`1`+ ches.data.wg$`0`
ches.data.wg$lrps.re <- rescale(ches.data.wg$lrps, to=c(0,1))

# merge with main dataset
parsys_ches <- left_join(parsys, ches.data.wg, by = c("country_name", "year")) 

ches_psdi <- ggplot(parsys_ches, aes(x = v2ps_democracy_adj, y =lrps.re))+
  geom_point()+
  geom_smooth(method = "lm", se=T)+
  xlab("Party-System Democracy Index")+
  ylab("Left-Right General Scale at Party-System Level (CHES, 1999-2019)") +
  stat_cor(label.x = 0, label.y = 0.90 )+
   labs(title = "A. Party Systems in Europe") +
  theme_light()+
  theme(legend.position = "bottom")

## Party System Closure
closure <- read_excel("casbertoa_Closure-2022.xlsx")
closure$country_name <- closure$country
parsys_closure <- left_join(parsys,closure, by = c("country_name", "year"))

psdi_closure <- ggplot(subset(parsys_closure,v2ps_democracy_adj>0),
                       aes(x = v2ps_democracy_adj, y = closure)) + 
  geom_point() +
  xlab("Party-System Democracy Index")+
  ylab("Party System Closure (Casal Bertoa and Eynedi, 2022)") +
  scale_x_continuous(breaks=seq(-1, 1, by = 0.25)) +
  labs(title= "B.") +
  stat_cor(label.x = 0, label.y = 50)+
  theme_light()+
  theme(legend.position = "bottom") 

psdi_closure <-psdi_closure + stat_smooth(method = "lm",
                                          formula = y ~ poly(x, 2))
### Parliamentary vs Presidential Systems ###
PPD <- read_excel("PPD.xls")

PPD$country_name <- PPD$countryname

PPD$pres.prl <- NA
PPD$pres.prl[PPD$legelec == 1 & PPD$preselec == 0] <- 0
PPD$pres.prl[PPD$legelec == 1 & PPD$preselec == 1] <- 1

compdata.cbb <- left_join(parsys, PPD,by = c("country_name", "year"))

ppd.plot <- ggplot(subset(compdata.cbb, !is.na(pres.prl) & v2ps_democracy_adj >0), 
                   aes(x = factor(pres.prl), y = v2ps_democracy_adj, 
                       fill = factor(pres.prl)))+
  geom_boxplot()+
  geom_signif(comparisons = list(c("0", "1")), 
              map_signif_level=TRUE)+
  geom_jitter(color="black", size=0.6, alpha=0.3) +
  xlab("Parliamentary & Presidential Systems (Cheibub, 2006)")+
  ylab("Party-System Democracy Index") +
  scale_y_continuous(breaks=seq(-1, 1, by = .25)) +
  labs(title = "C. ") +
  stat_cor(label.x = -0.7, label.y = -0.90 )+
  scale_fill_manual(values = wes_palette("BottleRocket2", n = 2),
                    labels = c("Parliamentary", "Presidential")) +
  labs(fill = "System: ") +
  theme_light()+
  theme(legend.position = "bottom")

(ches_psdi + psdi_closure)/
  ppd.plot+
  plot_layout(heights = c(2, 1))

## Figure A14. Cumulative Distribution Function ####
vparty <- vparty %>%
  filter(year >1969)%>%
  group_by(country_name) %>%
  mutate(Diff = lead(year) - year) %>%
  fill(Diff)

vparty$Diff <- ifelse(vparty$Diff>-1, vparty$Diff, NA)

ggplot(vparty, aes(x = Diff)) +
  stat_function(fun = pnorm, color = "black", size = 1) +
  theme_minimal() +
  labs(x = "Years Between Observations, V-Party",
       y = "Cumulative Probability")+
  xlim(1, 10)+
  scale_x_continuous(breaks=seq(0, 10, by = 1))

## Table A2. Full Results for Table 1 in Main Text ####
## Run line 354 on script "psdi_replication_ddl_and_empirics.R"
msummary(list_psdi_changes,estimate = "{estimate}{stars}", vcov = ~country_name)

## Table A3. Regime Change and PSDI, Global 1970-2019 ####
## Run line 492 on script "psdi_replication_ddl_and_empirics.R"
msummary(list_psdi_changes,estimate = "{estimate}{stars}", vcov = ~country_name)

## Table A4. Autocratization Onset around Election Time, 1970-2019 ####
psdi_ert_vdem_maddison_marquardt$reg_duration <- 
  psdi_ert_vdem_maddison_marquardt$reg_end_year-psdi_ert_vdem_maddison_marquardt$reg_start_year

psdi_ert_vdem_maddison_marquardt_a <- psdi_ert_vdem_maddison_marquardt%>%
  group_by(country_name, reg_duration)%>%
  mutate(reg_count = row_number())

## More Models
aut1_psdi <- glm(aut_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2)+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   (1|year), 
                 family = binomial(link = "probit"), psdi_ert_vdem_maddison_marquardt_a)

aut2_psdi <- glm(aut_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2)+
                   gdp_growth+gdppc_change+ln_gdppc+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   (1|year), 
                 family = binomial(link = "probit"), psdi_ert_vdem_maddison_marquardt_a)

aut3_psdi <- glm(aut_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2)+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   gdp_growth+gdppc_change+ln_gdppc+
                   lag.avg.dem.reg+lag.v2elaccept+
                   lag.v2xlg_legcon+
                   lag.v2x_cspart + z + (1|year), 
                 family = binomial(link = "probit"), psdi_ert_vdem_maddison_marquardt_a)

list.psdi_still <- list(aut1_psdi, aut2_psdi, aut3_psdi)
msummary(list.psdi_still,estimate = "{estimate}{stars}", vcov = "robust")

aut1_psdi <- glm(aut_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2)+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   (1|year), 
                 family = binomial(link = "probit"), psdi_ert_vdem_maddison_marquardt_a)

aut2_psdi <- glm(aut_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2)+
                   gdp_growth+gdppc_change+ln_gdppc+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   (1|year), 
                 family = binomial(link = "probit"), psdi_ert_vdem_maddison_marquardt_a)

aut3_psdi <- glm(aut_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2)+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   gdp_growth+gdppc_change+ln_gdppc+
                   lag.avg.dem.reg+lag.v2elaccept+
                   lag.v2xlg_legcon+
                   lag.v2x_cspart + z + (1|year), 
                 family = binomial(link = "probit"), psdi_ert_vdem_maddison_marquardt_a)

list.psdi_still <- list(aut1_psdi, aut2_psdi, aut3_psdi)
msummary(list.psdi_still,estimate = "{estimate}{stars}", vcov = "robust")

## Table A5. Democratization Onset around Election Time, 1970-2019 ####
dem1_psdi <- glm(dem_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2)+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   (1|year), 
                 family = binomial(link = "probit"), 
                 psdi_ert_vdem_maddison_marquardt_a)

dem2_psdi <- glm(dem_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2)+
                   gdp_growth+gdppc_change+ln_gdppc+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   (1|year), 
                 family = binomial(link = "probit"), 
                 psdi_ert_vdem_maddison_marquardt_a)

dem3_psdi <- glm(dem_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2)+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   gdp_growth+gdppc_change+ln_gdppc+
                   lag.avg.dem.reg+lag.v2elaccept+
                   lag.v2xlg_legcon+
                   lag.v2x_cspart + z + (1|year), 
                 family = binomial(link = "probit"), 
                 psdi_ert_vdem_maddison_marquardt_a)

list.psdi_demstill <- list(dem1_psdi, dem2_psdi, dem3_psdi)
msummary(list.psdi_demstill,estimate = "{estimate}{stars}", vcov = "robust")

## Table A6. Autocratization Onset predicted by Replacing PSDI with only Gvt. Democracy Index, 1970-2019 ####
aut1_psdi_gov <- glm(aut_onset_ele ~ v2ps_democracy_government+
                       reg_count+ I(reg_count^2)+I(reg_count^3)+
                       (1|year), 
                     family = binomial(link = "probit"), 
                     psdi_ert_vdem_maddison_marquardt_a)

aut2_psdi_gov <- glm(aut_onset_ele ~ v2ps_democracy_government+
                       gdp_growth+gdppc_change+ln_gdppc+
                       reg_count+ I(reg_count^2)+I(reg_count^3)+
                       (1|year), 
                     family = binomial(link = "probit"), 
                     psdi_ert_vdem_maddison_marquardt_a)

aut3_psdi_gov <- glm(aut_onset_ele ~ v2ps_democracy_government+
                       reg_count+ I(reg_count^2)+I(reg_count^3)+
                       gdp_growth+gdppc_change+ln_gdppc+
                       lag.avg.dem.reg+lag.v2elaccept+
                       lag.v2xlg_legcon+
                       lag.v2x_cspart + z + (1|year), 
                     family = binomial(link = "probit"), 
                     psdi_ert_vdem_maddison_marquardt_a)

list.psdi_gvtonly <- list(aut1_psdi_gov, aut2_psdi_gov, aut3_psdi_gov)
msummary(list.psdi_gvtonly,estimate = "{estimate}{stars}", vcov = "robust")

## Table A7. Democratization Onset predicted by Replacing PSDI with only Gvt. Democracy Index, 1970-2019 ####
dem1_psdi_gov <- glm(dem_onset_ele ~ v2ps_democracy_government+
                       reg_count+ I(reg_count^2)+I(reg_count^3)+
                       (1|year), 
                     family = binomial(link = "probit"), 
                     psdi_ert_vdem_maddison_marquardt_a)

dem2_psdi_gov <- glm(dem_onset_ele ~ v2ps_democracy_government+
                       gdp_growth+gdppc_change+ln_gdppc+
                       reg_count+ I(reg_count^2)+I(reg_count^3)+
                       (1|year), 
                     family = binomial(link = "probit"), 
                     psdi_ert_vdem_maddison_marquardt_a)

dem3_psdi_gov <- glm(dem_onset_ele ~ v2ps_democracy_government+
                       reg_count+ I(reg_count^2)+I(reg_count^3)+
                       gdp_growth+gdppc_change+ln_gdppc+
                       lag.avg.dem.reg+lag.v2elaccept+
                       lag.v2xlg_legcon+
                       lag.v2x_cspart + z + (1|year), 
                     family = binomial(link = "probit"), 
                     psdi_ert_vdem_maddison_marquardt_a)

list.psdi_gvtonly_dem <- list(dem1_psdi_gov, dem2_psdi_gov, dem3_psdi_gov)
msummary(list.psdi_gvtonly_dem,estimate = "{estimate}{stars}", vcov = "robust")

## Table A8. Democratic Tension and Regime Change, 1970-2019 ####

# creating a variable where positives is PSDI>EDI and negatives PSDI<EDI
vdem_multi$dem_tension1 <- vdem_multi$v2ps_democracy_2 - vdem_multi$v2x_polyarchy.x
## We have to account for regime change onsets with some room for years 
## before/after. This is because there could be cases where the PSDI
## observations are missing and replaced with the previous value BUT
## the EDI is changing constantly. So to limit possible drawbacks, we
## can allow for onset to exist within -1 to + 3 years, which would grant
## that if PSDI is *lower* than EDI at time t and and autocratization onset
## exists at time t+2 but no new data exist for PSDI (i.e., PSDI at t+2 == PSDI at t)
## then it is reasonable to expect that PSDI was going down and kept it lower
## than EDI. Although some errors might still exist, it seems to be a reasonable
## approach to assess whether the tension PSDI vs. EDI is a valuable one to 
## predict regime changes. If a party system is more/less democratic than the 
## regime, then it seems that the party systems can be driving the regime's 
## democratization/autocratization of the country. 


## Note: run script "psdi_replication_ddl_and_empirics lines 1 to 277 first
psdi_tension_aut <- glm(aut_onset_2yrs ~ dem_tension1+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        vdem_multi)

psdi_tension_aut2 <- glm(aut_onset_2yrs ~ dem_tension1+
                           gdp_growth+gdppc_change+ln_gdppc+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         vdem_multi)

psdi_tension_aut3 <- glm(aut_onset_2yrs ~ dem_tension1+
                           gdp_growth+gdppc_change+ln_gdppc+
                           lag.avg.dem.reg+lag.v2elaccept+lag.v2xlg_legcon+
                           lag.v2x_cspart+z+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         vdem_multi)

psdi_tension_dem <- glm(dem_onset_2yrs ~ dem_tension1+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        vdem_multi)

psdi_tension_dem2 <- glm(dem_onset_2yrs ~ dem_tension1+
                           gdp_growth+gdppc_change+ln_gdppc+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         vdem_multi)

psdi_tension_dem3 <- glm(dem_onset_2yrs ~dem_tension1+
                           gdp_growth+gdppc_change+ln_gdppc+
                           lag.avg.dem.reg+lag.v2elaccept+lag.v2xlg_legcon+
                           lag.v2x_cspart+z+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         vdem_multi)

list_psdi_tension <- list( psdi_tension_aut,psdi_tension_aut2,psdi_tension_aut3,
                           psdi_tension_dem,psdi_tension_dem2,psdi_tension_dem3)
msummary(list_psdi_tension,estimate = "{estimate}{stars}", vcov = "robust")

## Table A9. Main Results Dropping Cases Affected by the Collapse of Soviet Union ####
vdem_multi_noussr <- vdem_multi%>%
  filter(e_regionpol_6C != 1 & !between(year, 1989, 1999))

psdi_aut_nourrs <- glm(aut_ep_onset ~ v2ps_democracy_2+I(v2ps_democracy_2^2)+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         (1|year), 
                       family = binomial(link = "probit"), 
                       vdem_multi_noussr)

psdi_aut2_nourrs <- glm(aut_ep_onset ~ v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                          gdp_growth+gdppc_change+ln_gdppc+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        vdem_multi_noussr)

psdi_aut3_nourrs <- glm(aut_ep_onset ~ v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                          gdp_growth+gdppc_change+ln_gdppc+
                          lag.avg.dem.reg+lag.v2elaccept+lag.v2xlg_legcon+
                          lag.v2x_cspart+z+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        vdem_multi_noussr)

psdi_dem_nourrs <- glm(dem_ep_onset ~ v2ps_democracy_2+I(v2ps_democracy_2^2)+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         (1|year), 
                       family = binomial(link = "probit"), 
                       vdem_multi_noussr)

psdi_dem2_nourrs <- glm(dem_ep_onset ~ v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                          gdp_growth+gdppc_change+ln_gdppc+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        vdem_multi_noussr)

psdi_dem3_nourrs <- glm(dem_ep_onset ~ v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                          gdp_growth+gdppc_change+ln_gdppc+
                          lag.avg.dem.reg+lag.v2elaccept+lag.v2xlg_legcon+
                          lag.v2x_cspart+z+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        vdem_multi_noussr)

list_nourrs <- list(psdi_aut_nourrs, psdi_aut2_nourrs, psdi_aut3_nourrs,
                    psdi_dem_nourrs, psdi_dem2_nourrs, psdi_dem3_nourrs)
msummary(list_nourrs,estimate = "{estimate}{stars}", vcov = "robust")

## Table A10. PSDI and Regime Change Using Multinomial Logistic Regressions ####
vdem_multi$regime_change <- "No Regime Change"
vdem_multi$regime_change[vdem_multi$aut_ep_onset == 1] <- "Autocratization"
vdem_multi$regime_change[vdem_multi$dem_ep_onset == 1] <- "Democratization"

psdi_multinom <- multinom(regime_change ~ v2ps_democracy_2+I(v2ps_democracy_2^2)+
                            reg_count+ I(reg_count^2)+I(reg_count^3)+
                            (1|year), 
                          subset(vdem_multi, v2ps_democracy_2>0), Hess = TRUE)

psdi_multinom2 <- multinom(regime_change ~ v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                             gdp_growth+gdppc_change+ln_gdppc+
                             reg_count+ I(reg_count^2)+I(reg_count^3)+
                             (1|year),
                           subset(vdem_multi, v2ps_democracy_2>0), model = T)

psdi_multinom3 <- multinom(regime_change ~ v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                             gdp_growth+gdppc_change+ln_gdppc+
                             lag.avg.dem.reg+lag.v2elaccept+lag.v2xlg_legcon+
                             lag.v2x_cspart+z+
                             reg_count+ I(reg_count^2)+I(reg_count^3)+
                             (1|year), 
                           subset(vdem_multi, v2ps_democracy_2>0), model = T)


list_psdi_multinom <- list( psdi_multinom,psdi_multinom2,psdi_multinom3)
msummary(list_psdi_multinom,estimate = "{estimate}{stars}")
stargazer(list_psdi_multinom)

## Table A11. Cox Semi-parametric Analysis Replicating Main Results ####
vdem_multi_srv <- vdem_multi %>% 
  mutate(period = cumsum(aut_ep_onset == 1)) %>% ## which rows belong to one game period 
  group_by(period, country_name) %>% 
  mutate(aut_plus = case_when(aut_ep_onset == 1 ~ 0L , TRUE ~ row_number()),
         dem_plus = case_when(dem_ep_onset == 1 ~ 0L , TRUE ~ row_number())) %>%
  group_by(period,country_name) %>%  
  mutate(units = n() ) %>% ## how many rows per  game period
  mutate(aut_minus =  case_when(aut_ep_onset == 1 ~ 0L, TRUE ~  units - row_number() + 1L ),
         dem_minus =  case_when(dem_ep_onset == 1 ~ 0L, TRUE ~  units - row_number() + 1L )) %>%
  ungroup()


# check countries that had an autocratization onset
table(vdem_multi_srv$aut_ep_onset ==1, 
      vdem_multi_srv$country_name)

# Only subset to those countries that had an autocratization episode
vdem_multi_srv_aut <- vdem_multi_srv[vdem_multi_srv$country_name %in% 
                                       c("Afghanistan", "Algeria", "Argentina", "Armenia", "Azerbaijan",
                                         "Bahrain", "Bangladesh", "Belarus", "Benin", "Bolivia", "Botswana", "Brazil", "Bulgaria", 
                                         "Burkina Faso", "Burundi", "Cambodia", "Central African Republic",
                                         "Chad", "Chile", "Comoros", "Croatia", "Dominican Republic", "Ecuador", 
                                         "Egypt", "El Salvador", "Equatorial Guinea", "Estonia", "Eswatini", "Fiji", "Ghana",
                                         "Guatemala", "Guinea", "Guinea-Bissau", "Guyana", "Haiti", "Honduras", "Hong Kong",
                                         "Hungary", "India", "Indonesia", "Ivory Coast", "Kuwait",
                                         "Kyrgyzstan", "Laos", "Lesotho", "Liberia", "Libya", "Madagascar", "Malawi", "Maldives",
                                         "Mali", "Mauritania", "Mauritius", "Moldova", "Mongolia", "Nepal", "Nicaragua",
                                         "Niger", "Nigeria", "North Macedonia", "Pakistan", "Palestine/West Bank", 
                                         "Papua New Guinea", "Peru", "Philippines", "Poland", 
                                         "Republic of the Congo", "Russia", "Rwanda", "Serbia", "Seychelles",
                                         "Slovenia", "Solomon Islands", "South Korea", "Sri Lanka", "Sudan", "Suriname",
                                         "Tajikistan", "Tanzania", "Thailand", "The Gambia", "Tunisia", "Turkey", "Uganda",
                                         "Ukraine", "United States of America", "Uruguay", "Venezuela", "Yemen", "Zambia",
                                         "Zimbabwe"),]
vdem_multi_srv_aut

surv_fit_aut <- coxph(Surv(aut_minus, aut_ep) ~ 
                        v2ps_democracy_2+I(v2ps_democracy_2^2)+
                        reg_count+ I(reg_count^2)+I(reg_count^3)+
                        (1|year),
                      data = vdem_multi_srv_aut)

surv_fit_aut1 <- coxph(Surv(aut_minus, aut_ep) ~ 
                         v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                         gdp_growth+gdppc_change+ln_gdppc+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         (1|year),
                       data = vdem_multi_srv_aut)

surv_fit_aut2 <- coxph(Surv(aut_minus, aut_ep) ~ 
                         v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                         gdp_growth+gdppc_change+ln_gdppc+
                         lag.avg.dem.reg+lag.v2elaccept+lag.v2xlg_legcon+
                         lag.v2x_cspart+z+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         (1|year),
                       data = vdem_multi_srv_aut)

## Democratization
# check countries that had an autocratization onset
table(vdem_multi_srv$dem_ep_onset ==1, 
      vdem_multi_srv$country_name)

# Only subset to those countries that had an autocratization episode
vdem_multi_srv_dem <- vdem_multi_srv[vdem_multi_srv$country_name %in% 
                                       c("Afghanistan", "Albania", "Algeria", "Angola", "Argentina", "Armenia", "Azerbaijan",
                                         "Bahrain", "Bangladesh", "Belarus", "Benin", "Bhutan", "Bolivia", "Bosnia and Herzegovina",
                                         "Brazil", "Bulgaria", "Burkina Faso", "Burma/Myanmar", "Burundi", "Cambodia",
                                         "Cameroon", "Cape Verde", "Central African Republic", "Chad", 
                                         "Chile", "Colombia", "Comoros", "Costa Rica", "Croatia", "Cyprus",
                                         "Czechia", "Democratic Republic of the Congo", "Dominican Republic",
                                         "Ecuador", "El Salvador", "Equatorial Guinea", "Estonia", "Ethiopia",
                                         "Fiji", "Gabon", "Georgia", "German Democratic Republic", "Ghana", "Greece",
                                         "Guatemala", "Guinea", "Guinea-Bissau", "Guyana", "Haiti", "Honduras",
                                         "Hong Kong", "Hungary", "India", "Indonesia", "Iraq", "Ivory Coast",
                                         "Jamaica", "Jordan", "Kenya", "Kosovo", "Kuwait", "Kyrgyzstan", "Latvia",
                                         "Lebanon", "Lesotho", "Liberia", "Libya", "Madagascar","Malawi",  "Malaysia",
                                         "Maldives", "Mali", "Malta", "Mauritania", "Mexico", "Moldova", "Mongolia", "Montenegro",
                                         "Mozambique", "Namibia", "Nepal", "Nicaragua", "Niger","Nigeria", 
                                         "North Macedonia", "Oman", "Pakistan", "Palestine/West Bank", "Panama", "Papua New Guinea",
                                         "Paraguay", "Peru", "Philippines", "Poland", "Portugal", "Republic of the Congo",
                                         "Romania", "Russia", "Rwanda", "Sao Tome and Principe", "Senegal", "Serbia", "Seychelles",
                                         "Sierra Leone", "Slovakia", "Slovenia", "Solomon Islands",  "Somaliland",
                                         "South Africa", "South Korea", "Spain", "Sri Lanka", "Sudan",
                                         "Suriname", "Sweden", "Switzerland", "Taiwan", "Tanzania", "Thailand",
                                         "The Gambia", "Timor-Leste", "Togo", "Trinidad and Tobago", "Tunisia", "Turkey",
                                         "Uganda", "Ukraine","Uruguay", "Vanuatu", "Yemen", "Zambia", "Zanzibar", "Zimbabwe"),]

surv_fit_dem <- coxph(Surv(dem_minus, dem_ep) ~ 
                        v2ps_democracy_2+I(v2ps_democracy_2^2)+
                        reg_count+ I(reg_count^2)+I(reg_count^3)+
                        (1|year),
                      data = vdem_multi_srv_aut)

surv_fit_dem1 <- coxph(Surv(dem_minus, dem_ep) ~ 
                         v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                         gdp_growth+gdppc_change+ln_gdppc+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         (1|year),
                       data = vdem_multi_srv_aut)

surv_fit_dem2 <- coxph(Surv(dem_minus, dem_ep) ~ 
                         v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                         gdp_growth+gdppc_change+ln_gdppc+
                         lag.avg.dem.reg+lag.v2elaccept+lag.v2xlg_legcon+
                         lag.v2x_cspart+z+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         (1|year),
                       data = vdem_multi_srv_aut)

list_surv_all <- list(surv_fit_aut,surv_fit_aut1,surv_fit_aut2,
                      surv_fit_dem,surv_fit_dem1,surv_fit_dem2)
msummary(list_surv_all,estimate = "{estimate}{stars}", vcov = ~country_name)

## Figure A15. Cumulative Probabilities of Autocratization from Survival Model with PSDI ####
plot(survfit(surv_fit_aut2), #xlim=c(1,50), ylim=c(0.5,1),  
     xlab="Years Before Autocratization",
     ylab="Cumulative Probabilities of Autocratization",
     fun="event")

## Figure A16. Cumulative Probabilities of Democratization from Survival Model with PSDI ####
plot(survfit(surv_fit_dem2),  
     xlab="Years Before Democratization",
     ylab="Cumulative Probabilities of Democratization",
     fun="event")

## Table A12. Autocratization Onset and Party-System Institutionalization Measure, 1970-2019 ####

closure <- read_excel("casbertoa_Closure-2022.xlsx")
closure$country_name <- closure$country
psi_base <- read_csv("psi_base.csv")

psdi_psimeasures <- left_join(psdi_ert_vdem_maddison_marquardt_a,closure, by = c("country_name", "year"))%>%
  left_join(.,psi_base, by = c("country_name", "year"))

aut1_v2xps_party <- glm(aut_onset_ele ~ v2xps_party+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_psimeasures)

aut2_v2xps_party <- glm(aut_onset_ele ~ v2xps_party+
                          gdp_growth+gdppc_change+ln_gdppc+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_psimeasures)

aut3_v2xps_party <- glm(aut_onset_ele ~ v2xps_party+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          gdp_growth+gdppc_change+ln_gdppc+
                          lag.avg.dem.reg+lag.v2elaccept+
                          lag.v2xlg_legcon+
                          lag.v2x_cspart + z + (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_psimeasures)

aut1_closure <- glm(aut_onset_ele ~ closure+
                      reg_count+ I(reg_count^2)+I(reg_count^3)+
                      (1|year), 
                    family = binomial(link = "probit"), 
                    psdi_psimeasures)

aut2_closure <- glm(aut_onset_ele ~ closure+
                      gdp_growth+gdppc_change+ln_gdppc+
                      reg_count+ I(reg_count^2)+I(reg_count^3)+
                      (1|year), 
                    family = binomial(link = "probit"), 
                    psdi_psimeasures)

aut3_closure <- glm(aut_onset_ele ~ closure+
                      reg_count+ I(reg_count^2)+I(reg_count^3)+
                      gdp_growth+gdppc_change+ln_gdppc+
                      lag.avg.dem.reg+lag.v2elaccept+
                      lag.v2xlg_legcon+
                      (1|year), 
                    family = binomial(link = "probit"), 
                    psdi_psimeasures)

aut1_psi <- glm(aut_onset_ele ~ psi+
                  reg_count+ I(reg_count^2)+I(reg_count^3)+
                  (1|year), 
                family = binomial(link = "probit"), 
                psdi_psimeasures)

aut2_psi <- glm(aut_onset_ele ~ psi+
                  gdp_growth+gdppc_change+ln_gdppc+
                  reg_count+ I(reg_count^2)+I(reg_count^3)+
                  (1|year), 
                family = binomial(link = "probit"), 
                psdi_psimeasures)

aut3_psi <- glm(aut_onset_ele ~ psi+
                  reg_count+ I(reg_count^2)+I(reg_count^3)+
                  gdp_growth+gdppc_change+ln_gdppc+
                  lag.avg.dem.reg+lag.v2elaccept+
                  lag.v2xlg_legcon+
                  lag.v2x_cspart + z + (1|year), 
                family = binomial(link = "probit"), 
                psdi_psimeasures)

list.psi <- list(aut1_v2xps_party, aut2_v2xps_party, aut3_v2xps_party,
                 aut1_closure, aut2_closure, aut3_closure,
                 aut1_psi, aut2_psi, aut3_psi)
msummary(list.psi,estimate = "{estimate}{stars}", vcov = "robust")
#
## Table A13. Democratization Onset and Party-System Institutionalization Measure, 1970-2019 ####
dem1_v2xps_party <- glm(dem_onset_ele ~ v2xps_party+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_psimeasures)

dem2_v2xps_party <- glm(dem_onset_ele ~ v2xps_party+
                          gdp_growth+gdppc_change+ln_gdppc+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_psimeasures)

dem3_v2xps_party <- glm(dem_onset_ele ~ v2xps_party+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          gdp_growth+gdppc_change+ln_gdppc+
                          lag.avg.dem.reg+lag.v2elaccept+
                          lag.v2xlg_legcon+
                          lag.v2x_cspart + z + (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_psimeasures)

dem1_closure <- glm(dem_onset_ele ~ closure+
                      reg_count+ I(reg_count^2)+I(reg_count^3)+
                      (1|year), 
                    family = binomial(link = "probit"), 
                    psdi_psimeasures)

dem2_closure <- glm(dem_onset_ele ~ closure+
                      gdp_growth+gdppc_change+ln_gdppc+
                      reg_count+ I(reg_count^2)+I(reg_count^3)+
                      (1|year), 
                    family = binomial(link = "probit"), 
                    psdi_psimeasures)

dem3_closure <- glm(dem_onset_ele ~ closure+
                      reg_count+ I(reg_count^2)+I(reg_count^3)+
                      gdp_growth+gdppc_change+ln_gdppc+
                      lag.avg.dem.reg+lag.v2elaccept+
                      lag.v2xlg_legcon+
                      lag.v2x_cspart + z + (1|year), 
                    family = binomial(link = "probit"), 
                    psdi_psimeasures)

dem1_psi <- glm(dem_onset_ele ~ psi+
                  reg_count+ I(reg_count^2)+I(reg_count^3)+
                  (1|year), 
                family = binomial(link = "probit"), 
                psdi_psimeasures)

dem2_psi <- glm(dem_onset_ele ~ psi+
                  gdp_growth+gdppc_change+ln_gdppc+
                  reg_count+ I(reg_count^2)+I(reg_count^3)+
                  (1|year), 
                family = binomial(link = "probit"), 
                psdi_psimeasures)

dem3_psi <- glm(dem_onset_ele ~ psi+
                  reg_count+ I(reg_count^2)+I(reg_count^3)+
                  gdp_growth+gdppc_change+ln_gdppc+
                  lag.avg.dem.reg+lag.v2elaccept+
                  lag.v2xlg_legcon+
                  lag.v2x_cspart + z + (1|year), 
                family = binomial(link = "probit"), 
                psdi_psimeasures)

list.psi_dem <- list(dem1_v2xps_party, dem2_v2xps_party, dem3_v2xps_party,
                     dem1_closure, dem2_closure, dem3_closure,
                     dem1_psi, dem2_psi, dem3_psi)
msummary(list.psi_dem,estimate = "{estimate}{stars}", vcov = "robust")

## Table A14-A16. PSDI Changes over three election years and influence on regime change, 1970-2019 ####
## PSDI CHANGE 1-ELE -->AUTOCRATIZATION ###
aut1_psdi_change <- glm(aut_onset_ele ~ psdi_change1+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_ert_vdem_maddison_marquardt_a)

aut2_psdi_change <- glm(aut_onset_ele ~ psdi_change1+
                          gdp_growth+ln_gdppc+gdppc_change+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_ert_vdem_maddison_marquardt_a)

aut3_psdi_change <- glm(aut_onset_ele ~ psdi_change1+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          gdp_growth+ln_gdppc+gdppc_change+
                          lag.avg.dem.reg+lag.v2elaccept+
                          lag.v2xlg_legcon+
                          lag.v2x_cspart + z + (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_ert_vdem_maddison_marquardt_a)
## PSDI CHANGE 2-ELE -->AUTOCRATIZATION ###
aut1_psdi_change2 <- glm(aut_onset_ele ~ psdi_change2+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)

aut2_psdi_change2 <- glm(aut_onset_ele ~ psdi_change2+
                           gdp_growth+ln_gdppc+gdppc_change+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)

aut3_psdi_change2 <- glm(aut_onset_ele ~ psdi_change2+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           gdp_growth+ln_gdppc+gdppc_change+
                           lag.avg.dem.reg+lag.v2elaccept+
                           lag.v2xlg_legcon+
                           lag.v2x_cspart + z + (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)

## PSDI CHANGE 3-ELE -->AUTOCRATIZATION ###
aut1_psdi_change3 <- glm(aut_onset_ele ~ psdi_change3+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)

aut2_psdi_change3 <- glm(aut_onset_ele ~ psdi_change3+
                           gdp_growth+ln_gdppc+gdppc_change+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)

aut3_psdi_change3 <- glm(aut_onset_ele ~ psdi_change3+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           gdp_growth+ln_gdppc+gdppc_change+
                           lag.avg.dem.reg+lag.v2elaccept+
                           lag.v2xlg_legcon+
                           lag.v2x_cspart + z + (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)

list_aut_psdi_change1 <- list(aut1_psdi_change,aut2_psdi_change,aut3_psdi_change,
                              aut1_psdi_change2,aut2_psdi_change2,aut3_psdi_change2,
                              aut1_psdi_change3,aut2_psdi_change3,aut3_psdi_change3)
msummary(list_aut_psdi_change1,estimate = "{estimate}{stars}", vcov =  ~country_name)
#
## PSDI CHANGE 1-ELE -->DEMOCRATIZATION ###
dem1_psdi_change <- glm(dem_onset_ele ~ psdi_change1+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_ert_vdem_maddison_marquardt_a)

dem2_psdi_change <- glm(dem_onset_ele ~ psdi_change1+
                          gdp_growth+ln_gdppc+gdppc_change+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_ert_vdem_maddison_marquardt_a)

dem3_psdi_change <- glm(dem_onset_ele ~ psdi_change1+
                          reg_count+ I(reg_count^2)+I(reg_count^3)+
                          gdp_growth+ln_gdppc+gdppc_change+
                          lag.avg.dem.reg+lag.v2elaccept+
                          lag.v2xlg_legcon+
                          lag.v2x_cspart + z + 
                          (1|year), 
                        family = binomial(link = "probit"), 
                        psdi_ert_vdem_maddison_marquardt_a)
## PSDI CHANGE 2-ELE -->DEMOCRATIZATION ###
dem1_psdi_change2 <- glm(dem_onset_ele ~ psdi_change2+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)

dem2_psdi_change2 <- glm(dem_onset_ele ~ psdi_change2+
                           gdp_growth+ln_gdppc+gdppc_change+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)

dem3_psdi_change2 <- glm(dem_onset_ele ~ psdi_change2+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           gdp_growth+ln_gdppc+gdppc_change+
                           lag.avg.dem.reg+lag.v2elaccept+
                           lag.v2xlg_legcon+
                           lag.v2x_cspart + z + 
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)

## PSDI CHANGE 3-ELE -->DEMOCRATIZATION ###
dem1_psdi_change3 <- glm(dem_onset_ele ~ psdi_change3+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a,
                         na.action = na.exclude)

dem2_psdi_change3 <- glm(dem_onset_ele ~ psdi_change3+
                           gdp_growth+ln_gdppc+gdppc_change+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a,
                         na.action = na.omit)

dem3_psdi_change3 <- glm(dem_onset_ele ~ psdi_change3+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           gdp_growth+ln_gdppc+gdppc_change+
                           lag.avg.dem.reg+
                           lag.v2xlg_legcon+ lag.v2elaccept+
                           lag.v2x_cspart  + z+
                           (1|year),
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a,
                         na.action = na.omit)

list_dem_psdi_change1 <- list(dem1_psdi_change,dem2_psdi_change,dem3_psdi_change,
                              dem1_psdi_change2,dem2_psdi_change2,dem3_psdi_change2,
                              dem1_psdi_change3,dem2_psdi_change3,dem3_psdi_change3)
msummary(list_dem_psdi_change1,estimate = "{estimate}{stars}", vcov =  ~country_name)
#
## Table A17. % change in Democratic levels only gvt and autocratization, 1970-2019 ####
aut1_psdi_gvt_change <- glm(aut_onset_ele ~ psdi_gov_change+
                              reg_count+ I(reg_count^2)+I(reg_count^3)+
                              (1|year), 
                            family = binomial(link = "probit"), 
                            psdi_ert_vdem_maddison_marquardt_a)

aut2_psdi_gvt_change <- glm(aut_onset_ele ~ psdi_gov_change+
                              gdp_growth+ln_gdppc+gdppc_change+
                              reg_count+ I(reg_count^2)+I(reg_count^3)+
                              (1|year), 
                            family = binomial(link = "probit"), 
                            psdi_ert_vdem_maddison_marquardt_a)

aut3_psdi_gvt_change <- glm(aut_onset_ele ~ psdi_gov_change+
                              reg_count+ I(reg_count^2)+I(reg_count^3)+
                              gdp_growth+ln_gdppc+gdppc_change+
                              lag.avg.dem.reg+lag.v2elaccept+
                              lag.v2xlg_legcon+
                              lag.v2x_cspart + z + (1|year), 
                            family = binomial(link = "probit"), 
                            psdi_ert_vdem_maddison_marquardt_a)
list_aut_gvt_change <- list(aut1_psdi_gvt_change,aut2_psdi_gvt_change,aut3_psdi_gvt_change)
msummary(list_aut_gvt_change,estimate = "{estimate}{stars}", vcov =  ~country_name)

## Table A18. % Change in party systems’ indexes and autocratization ####
psdi_psimeasures <- psdi_psimeasures%>%
  group_by(country_name)%>%
  mutate(v2xps_party_change_a = ((v2xps_party / lag(v2xps_party) - 1) * 100),
         closure_change = ((closure / lag(closure) - 1) * 100),
         psi_change = ((psi / lag(psi) - 1) * 100)
  )

# Autocratization Onset BHS 2017 CHANGE ###
aut1_v2xps_party_change_change <- glm(aut_onset_ele ~ v2xps_party_change+
                                        reg_count+ I(reg_count^2)+I(reg_count^3)+
                                        (1|year), 
                                      family = binomial(link = "probit"), 
                                      psdi_psimeasures)

aut2_v2xps_party_change_change <- glm(aut_onset_ele ~ v2xps_party_change+
                                        gdp_growth+ln_gdppc+gdppc_change+
                                        reg_count+ I(reg_count^2)+I(reg_count^3)+
                                        (1|year), 
                                      family = binomial(link = "probit"), 
                                      psdi_psimeasures)

aut3_v2xps_party_change_change <- glm(aut_onset_ele ~ v2xps_party_change+
                                        reg_count+ I(reg_count^2)+I(reg_count^3)+
                                        gdp_growth+ln_gdppc+gdppc_change+
                                        lag.avg.dem.reg+lag.v2elaccept+
                                        lag.v2xlg_legcon+
                                        lag.v2x_cspart + z + (1|year), 
                                      family = binomial(link = "probit"), 
                                      psdi_psimeasures)
## Autocratization Onset CSE 2016 CHANGE ###
aut1_closure_change <- glm(aut_onset_ele ~ closure_change+
                             reg_count+ I(reg_count^2)+I(reg_count^3)+
                             (1|year), 
                           family = binomial(link = "probit"), 
                           psdi_psimeasures)

aut2_closure_change <- glm(aut_onset_ele ~ closure_change+
                             gdp_growth+ln_gdppc+gdppc_change+
                             reg_count+ I(reg_count^2)+I(reg_count^3)+
                             (1|year), 
                           family = binomial(link = "probit"), 
                           psdi_psimeasures)

aut3_closure_change <- glm(aut_onset_ele ~ closure_change+
                             reg_count+ I(reg_count^2)+I(reg_count^3)+
                             gdp_growth+ln_gdppc+gdppc_change+
                             lag.avg.dem.reg+lag.v2elaccept+
                             lag.v2xlg_legcon+
                             lag.v2x_cspart + z + (1|year), 
                           family = binomial(link = "probit"), 
                           psdi_psimeasures)
### Autocratization Onset Kim 2023 CHANGE ###
aut1_psi_change <- glm(aut_onset_ele ~ psi_change+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         (1|year), 
                       family = binomial(link = "probit"), 
                       psdi_psimeasures)

aut2_psi_change <- glm(aut_onset_ele ~ psi_change+
                         gdp_growth+ln_gdppc+gdppc_change+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         (1|year), 
                       family = binomial(link = "probit"), 
                       psdi_psimeasures)

aut3_psi_change <- glm(aut_onset_ele ~ psi_change+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         gdp_growth+ln_gdppc+gdppc_change+
                         lag.avg.dem.reg+lag.v2elaccept+
                         lag.v2xlg_legcon+
                         lag.v2x_cspart + z + (1|year), 
                       family = binomial(link = "probit"), 
                       psdi_psimeasures)

list.auto_ps <- list(
  aut1_v2xps_party_change_change,aut2_v2xps_party_change_change,aut3_v2xps_party_change_change,
  aut1_closure_change,aut2_closure_change,aut3_closure_change,
  aut1_psi_change, aut2_psi_change, aut3_psi_change
)
msummary(list.auto_ps,estimate = "{estimate}{stars}", vcov = ~country_name)

## Table A19. Government Dem. % and Democratization ####
dem1_psdi_gvt_change <- glm(dem_onset_ele ~ psdi_gov_change+
                              reg_count+ I(reg_count^2)+I(reg_count^3)+
                              (1|year), 
                            family = binomial(link = "probit"), 
                            psdi_ert_vdem_maddison_marquardt_a)

dem2_psdi_gvt_change <- glm(dem_onset_ele ~ psdi_gov_change+
                              gdp_growth+ln_gdppc+gdppc_change+
                              reg_count+ I(reg_count^2)+I(reg_count^3)+
                              (1|year), 
                            family = binomial(link = "probit"), 
                            psdi_ert_vdem_maddison_marquardt_a)

dem3_psdi_gvt_change <- glm(dem_onset_ele ~ psdi_gov_change+
                              reg_count+ I(reg_count^2)+I(reg_count^3)+
                              gdp_growth+ln_gdppc+gdppc_change+
                              lag.avg.dem.reg+lag.v2elaccept+
                              lag.v2xlg_legcon+
                              lag.v2x_cspart + z + (1|year), 
                            family = binomial(link = "probit"), 
                            psdi_ert_vdem_maddison_marquardt_a)
list_dem_gvt_change <- list(dem1_psdi_gvt_change,dem2_psdi_gvt_change,dem3_psdi_gvt_change)
msummary(list_dem_gvt_change,estimate = "{estimate}{stars}", vcov =  ~country_name)

## Table A20. Democratization and % change in party systems measures #### 

# Democratization Onset BHS 2017 CHANGE ###
dem1_v2xps_party_change_change <- glm(dem_onset_ele ~ v2xps_party_change+
                                        reg_count+ I(reg_count^2)+I(reg_count^3)+
                                        (1|year), 
                                      family = binomial(link = "probit"), 
                                      psdi_psimeasures)

dem2_v2xps_party_change_change <- glm(dem_onset_ele ~ v2xps_party_change+
                                        gdp_growth+ln_gdppc+gdppc_change+
                                        reg_count+ I(reg_count^2)+I(reg_count^3)+
                                        (1|year), 
                                      family = binomial(link = "probit"), 
                                      psdi_psimeasures)

dem3_v2xps_party_change_change <- glm(dem_onset_ele ~ v2xps_party_change+
                                        reg_count+ I(reg_count^2)+I(reg_count^3)+
                                        gdp_growth+ln_gdppc+gdppc_change+
                                        lag.avg.dem.reg+lag.v2elaccept+
                                        lag.v2xlg_legcon+
                                        lag.v2x_cspart + z + (1|year), 
                                      family = binomial(link = "probit"), 
                                      psdi_psimeasures)
## Democratization Onset CSE 2016 CHANGE ###
dem1_closure_change <- glm(dem_onset_ele ~ closure_change+
                             reg_count+ I(reg_count^2)+I(reg_count^3)+
                             (1|year), 
                           family = binomial(link = "probit"), 
                           psdi_psimeasures)

dem2_closure_change <- glm(dem_onset_ele ~ closure_change+
                             gdp_growth+ln_gdppc+gdppc_change+
                             reg_count+ I(reg_count^2)+I(reg_count^3)+
                             (1|year), 
                           family = binomial(link = "probit"), 
                           psdi_psimeasures)

dem3_closure_change <- glm(dem_onset_ele ~ closure_change+
                             reg_count+ I(reg_count^2)+I(reg_count^3)+
                             gdp_growth+ln_gdppc+gdppc_change+
                             lag.avg.dem.reg+lag.v2elaccept+
                             lag.v2xlg_legcon+
                             lag.v2x_cspart + z + (1|year), 
                           family = binomial(link = "probit"), 
                           psdi_psimeasures)
### Democratization Onset Kim 2023 CHANGE ###
dem1_psi_change <- glm(dem_onset_ele ~ psi_change+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         (1|year), 
                       family = binomial(link = "probit"), 
                       psdi_psimeasures)

dem2_psi_change <- glm(dem_onset_ele ~ psi_change+
                         gdp_growth+ln_gdppc+gdppc_change+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         (1|year), 
                       family = binomial(link = "probit"), 
                       psdi_psimeasures)

dem3_psi_change <- glm(dem_onset_ele ~ psi_change+
                         reg_count+ I(reg_count^2)+I(reg_count^3)+
                         gdp_growth+ln_gdppc+gdppc_change+
                         lag.avg.dem.reg+lag.v2elaccept+
                         lag.v2xlg_legcon+
                         lag.v2x_cspart + z + (1|year), 
                       family = binomial(link = "probit"), 
                       psdi_psimeasures)

list.auto_ps <- list(
  dem1_v2xps_party_change_change,dem2_v2xps_party_change_change,dem3_v2xps_party_change_change,
  dem1_closure_change,dem2_closure_change,dem3_closure_change,
  dem1_psi_change, dem2_psi_change, dem3_psi_change
)
msummary(list.auto_ps,estimate = "{estimate}{stars}", vcov = ~country_name)
